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Abstract. A stochastic algorithm for simulation of fluctuation-induced kinetics 
of H 2 formation on grain surfaces is suggested as a generalization of the 
technique developed in our recent studies [14] where this method was developed to 
describe the annihilation of spatially separate electrons and holes in a disordered 
semiconductor. The stochastic model is based on the spatially inhomogeneous, 
nonlinear integro-differential Smoluchowski equations with random source term. 
In this paper we derive the general system of Smoluchowski type equations for 
the formation of H 2 from two hydrogen atoms on the surface of interstellar dust 
grains with physisorption and chemisorption sites. We focus in this study on 
the spatial distribution, and numerically investigate the segregation in the case 
of a source with a continuous generation in time and randomly distributed in 
space. The stochastic particle method presented is based on a probabilistic 
interpretation of the underlying process as a stochastic Markov process of 
interacting particle system in discrete but randomly progressed time instances. 
The segregation is analyzed through the correlation analysis of the vector random 
field of concentrations which appears to be isotropic in space and stationary in 
time. 
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1. Introduction 

In the literature, there exist several theoretical models and methods to study the 
surface chemistry 

that occurs on interstellar dust grains under different astrophysical conditions. 
We mention here the following approaches: the conventional rate equation method 
which can be implemented using the master equation, or a Monte Carlo technique 
which are able to simulate the diffusive reactions on surfaces. The kinetic Monte 
Carlo technique has been applied in all branches of astrochemistry studies, see, e.g. 
05 0; 0) 0 and seems to be the most popular in practical calculations. The 
moment equation approach presents a nice analytical approach, however it is much 
weakened by the closer assumptions mu. 

The fluctuation impact on the kinetics in the interstellar chemistry has not 
been studied thoroughly before, in spite the fact that in chemical kinetics, the 
segregation has been studied by many authors, see, for instance [$j where the 
segregation phenomenon has been described from a general point of view in 
the case of a source, random in space and continuous in time. Note that the 
Monte Carlo methods used in studies mentioned above are of different type: they 
do not deal with the spatial fluctuations of the source, actually, they suggest 
a randomized solution of deterministic rate equations which do not govern the 
segregation phenomenon. 

Let explain this by a simple example by considering a reaction of two types 
of particles, A and B, leading to a product P: A + B —>■ P. The simplest kinetic 
approach to this reaction first considered by Smoluchowski [XBJ is based on the rate 
equations 

dnA/dt — —g dns/dt — —gn A (t)nB{t), ua( 0) = do, tib(O) = b 0 

where g is the reaction rate. For diffusion-controlled reactions in R 3 , Smoluchowski 
obtained g = AnDro, where ro is the particle radius and D is the relative diffusion 
coefficient. This equation can be easily solved explicitly: riA{t) = «o/(l + oo gt) for 
a 0 = b 0 , and n A (t) = a 0 fo/{b 0 exp(f 0 gt) - a 0 \ for a 0 < b 0 , where f 0 = b 0 - a 0 . 

Such description of chemical reactions implies conditions such that the rate at 
which the reactants approach each other (the diffusion rate) is much larger than the 
rate at which they react chemically. Hence, the basic assumption underlying this 
reaction model is a homogeneous spatial distribution of particles during the reaction 
at any time instant, i. e., the components A and B should be always perfectly 
mixed. Independent of the dimension d, this assumption results in solutions with 
a long-time asymptotics ~ exp(— fogt) if a 0 b 0 , and r\_/ 1 /gt if a 0 = b 0 . 

In this approach, the diffusion is treated macroscopically, ignoring density 
fluctuations. In the absence of nonlinear interactions such as chemical reactions, 
the macroscopic diffusion equations govern uniform concentration distributions. 
Thermal fluctuations, initial density inhomogeneities, and the randomness of 
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reaction events lead to non-uniform concentration fields and changes the solution 
structure drastically, in particular, the time dependence of the mean solution for 
asymptotically long times. 

Fluctuations are responsible for the spatial correlations, and in particular, 
they may lead to segregation, i.e., the formation of spatially separated clusters 
composed entirely of particles of either type A or B. In this paper we of appearing 
segregation in terms of the input parameters. 

It was first shown by Ovchinnikov and Zeldovich m that in the fluctuation- 
induced reactions, the long-time asymptotics is ~ f~ 3//4 if a 0 = b 0 . Generally, for 
ao = the asymptotics ~ f _d//4 is valid for any dimensionality d < 4, while for 
d > 4, ~ f _1 . This law was obtained by several authors using different arguments 
(see, e.g., Refs. mi, n, 0, and mi)- If the densities of the components are not 
equal, the asymptotics is different. For instance, if «o < n a ~ aoexp(— y/t) for 
d — 1, ua ~ clq exp {—t/ log(t)} for d — 2, while for d > 3, the asymptotic law again 
coincides with the homogeneous case: ua ~ aoexp(— t) [lj. 

This difference between the homogeneous and fluctuation-limited kinetics 
holds regardless of the type of reaction between the particles. It applies, for 
instance, to the Smoluchowski coagulation equation. Coagulation, or coalescence, 
is a process by which two particles collide and adhere, or coagulate. There are many 
different mechanisms that bring two particles to each other: molecular diffusion, 
gravitational sedimentation, free molecule collisions, turbulent motion of the host 
gas, acoustic waves, density, concentration and temperature gradients, electric 
charges, etc. (see, e.g., Refs. |19j and ED- Let us consider the Smoluchowski 
equation governing the coagulation of particles of different type having different 
masses mp (e.g., see BE|: 


dn(m, t) 
dt 


1 

2 


m i 


R(u, m — u)n(m — u, t)n(u , t)du 


— n(m,t ) J ... j R(u,m)n(u,t)du. (1) 
o o 

Here, rn t is the mass of the i-th component in a particle, and m is a vector of 
compositions (mi,... ,m s ), where s is the total number of components; n(m,t)dm 
is the number of particles having mass of component i in the range [rn t , m, + dmf\ 
at time f, and R(u,m) = R(in,u) is the binary coagulation coefficient. 

The numerical solution of the inhomogeneous Smoluchowski equations is a 
highly challenging problem even for only a few particle types. We deal in this 
paper with three particle types (physisorbed H (Hp, chemisorbed H (He), and 
molecular hydrogen H 2 ). 

The main difficulties in the problem we solve, can be formulated as follows: 

(1) The major difficulty arises from the inhomogeneity in space. 

(2) The second difficulty of our problem is caused by the low and singular particle 
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densities. 

(3) Third, we are interested in the particle kinetics for very long times, to reach 
and study a quasi-stationary regime. 

(4) Finally, we deal with stochastic source term for phisisorbed H, and we have to 
take the average over a reach ensemble of source realizations. In addition, other 
parameters of the governing equations may fluctuate randomly and have a large 
impact on the kinetics of the processes studied. 

Conventional numerical methods are not applicable for handling problems of 
this kind, and we therefore use the Monte Carlo methods we have developed in our 
recent papers [2] and [45] which were based on the approach previously suggested 
in our publications Refs. [13], |7j, and ra. 

The main idea of Monte Carlo methods for solving the spatially homogeneous 
Smoluchowski equation lies in the probabilistic interpretation of the evolution of 
the interacting particles as a Markov chain (see, e. g., Ref. HD- In Refs. [6]| and ta. 
we have applied the Monte Carlo technique to the inhomogeneous Smoluchowski 
equation. In this paper, we consider the general inhomogeneous case with diffusion. 
Note that the direct Monte Carlo simulation of the particle interactions and 
diffusion jumps on a grid is computationally expensive, because one has to consider 
a huge number of jumps per one particle interaction. In [14] we suggest a new Monte 
Carlo method for this case, introducing ’’long diffusion jumps” which accelerate 
the simulation process significantly. However this approach being very efficient for 
calculating most important statistical characteristics of the solution, is not well 
adjusted when there is a need to analyze the spatial distribution, in particular, 
the segregation process. Therefore, we use a detailed Monte Carlo scheme where 
in each random time step, all processes including the diffusion jump are directly 
simulated. 

2. The inhomogeneous nonlinear Smoluchowski equation based model 

We assume that hydrogen atoms absorbed on a grain site can move to another site 
via tunneling or thermal diffusion. We define Hp , He and H- 2 as the physisorbed 
H, chemisorbed H and molecular hydrogen concentrations, respectively. We may 
imagine that the granular surfaces to be square lattices with four nearest neighbor 
sites, as on fcc[100] plane. But we would prefer to use a continuous description, both 
in space and time. This approach dramatically reduces the memory and computer 
time needed, and makes it possible to simulate the H 2 formation of very large grain 
surfaces. In the case of low temperatures, when the diffusion is negligible, we work 
on a grid. 

Different energy barriers occur between different pairs of sites which can be 
sites in which H atoms are weakly bound, due to physisorption (weak Van der 
Waals interaction), or strongly covalent bound, due to chemisorption. We assume 
the number of physisorbed and chemisorbed sites on a grain to be identical. 
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We denote by Dp the diffusion coefficient for the physisorption, and by Dc the 
chemisorption diffusion coefficient. The recombination of spatially separate sites 
distributed in a grain G due to tunneling from a physisorption to physisorption site 
is defined by the rate app(r), from a chemisorption to physisorption site it is defined 
by acp(j). Analogously, the tunneling from a chemisorption to chemisorption 
site is defined by the rate acc( r ), from a physisorption site to chemisorption 
site is defined by the rate apc(r). These tunneling coefficients have the form 
Qjj(r) = a^exp(—\x\/aij) where \x\ is the distance between the two interacting 
sites. Here is the characteristic distance of tunneling, and (i,j=P,C) are the 
frequencies of the relevant events which are described in details in [5], see also [2]. 

The spatially varying desorption rates are designated by Wh p and Wh c - The 
source F( r) is the accretion rate. In this study we assume that F( r) is a spatial 
random field, and in simplest case we use the Poissonian distribution. 

Altogether, the concentrations Hp(r,t ) and Hp(r,t ) are governed by the 
following system of two coupled nonlinear Smoluchowski equation 




Without loss of generality, we assume that the total initial concentrations of 
Hp and of Hq are equal. At the initial time t — 0, the concentrations Hp and He 
are zero. 

Since the source F is a random field, so are the concentrations Hp and Hq- 
Therefore, the experimentally measured is the mean concentration of the molecules 
H 2 . The kinetics of H 2 concentration reads 


dH 



We are interested both in kinetics and the quasi-steady state which is reached 
in which the total mean surface population of H atoms fluctuates around a constant 
value. After the steady-state has been reached the H 2 formation efficiency is then 
defined by the ratio of the means: 



(5) 


where the angle brackets stand for the ensemble average generated by the random 
source F. 
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3. Monte Carlo Algorithms 

The equations (§-!§ have the structure of inhomogeneous Smoluchowski 
coagulation equations, as mentioned in the introduction. The Smoluchowski 
equations can be interpreted probabilistically as an equation generated by Markov 
chains describing the evolution of pairwise interacting particle system. In Ref. 0 ). 
we developed a Monte Carlo algorithm for inhomogeneous Smoluchowski equation, 
which we adapt in 0 - Here it will be used to solve Eqs. (j2]{3]) for the two- 
dimensional case d — 2 , with the focus on the segregation problem. 


3.1. Recombinations by tunneling, in the absence of diffusion 

For simplicity, let us first consider the algorithm for the case when there is no 
diffusion. The direct simulation assumes that the process of interactions is pairwise 
and Markovian, i.e., having made a time step, the next time step is simulated 
independently. The interacting pair is sampled according to the kernel of the 
equation B = B 0 exp(— |x|/a) as described below. 

Let us assume our phase space is continuous. The process is simulated in a 
square box with size L x L and periodic boundary conditions. 

In the first step, we sample between the following possible events: (i) formation 
of H 2 by tunneling from Hp to Hp, (ii) formation of H 2 by tunneling from Hp to 
He, (Hi) formation of H 2 by tunneling from He to He, and ( iv ) formation of 
Ho by tunneling from He to Hp. The processes are described by the relevant 
probabilities ay,(r) = a- exp(— \x\/aif). Thus, the simulation algorithm in the 
absence of diffusion can be described as follows. 


(i) Put t = 0, and sample n = n 0 atoms Hp and p = n 0 atoms H c at random, 
independently and uniformly distributed in the box. 

(ii) Sample one of the possible events: (i) Hp +Hp , ( ii ) Hp +He , (in) He +He 
and (iv) He +Hp . To do this, first calculate the majorant frequencies for the 
four events: 


n(n - 1 ) 0 r 

Ai =--- a PP exp {- 


' nn.min 


a P p 


h 


X 2 = npot PC exp {- 
A 3 = npcf CP exp {- 


1 np,min 

a P c 

r np,min 


dCP 


x _p(p~ !)„ 0 

A4 —-—- oi, 


cc 


exp {- 


' pp,mm 

&CC 


h 


( 6 ) 

(7) 

( 8 ) 

(9) 


where r nP)m i n is the minimal of all possible distances between n Hp atoms and p 
He atoms in the box, and analogously for other cases. From these frequencies, 
calculate the probabilities p 1 , p 2 , pn , p 4 of the events (i), (ii), (Hi) and (iv), 
respectively: pi = Ai/A, p 2 = A 2 /A, p% = A 3 /A and p 4 = 1 — p\— p 2 — p 3 , where 
A = Ai + A 2 + A 3 + A 4 . 
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(iii) From the probabilities pi,P 2 ,P 3 ,P±, sample the event k = i,ii,iii,iv, calculate 
the time increment as A t = — \og{rand) / X, and calculate t := t + At. 

(iv) For the sampled event k{— i,ii,iii,iv), choose uniformly the relevant 
interacting pair, and check if the interaction takes place. For instance, if 
k = ii, i. e., the sampled event happens to be recombination Hp +Hc , 
calculate P np = exp { ~ r+r ’" p ' min }. If rand < P np , then the event occurs. Hence, 
recalculate n := n — 1 and p p — 1, and go back to step 2. Otherwise, 
nothing happens, the probabilities Pi,P 2 ,P 3 ,Pa remain the same, and so return 
to step 3. Here r is the distance between the sampled interacting pair. For the 
case k — i, put n := n — 2. If the concentrations n and p are calculated at some 
prescribed time instances t m , m — 1,... M, just score the values n(t m ), p(t m ), 
m = 1,... M. To calculate the concentration of H 2 , say, at a time t from 
the interval t G [UU 2 ], count M, the number of all recombinations which have 
occurred during this time interval, and take the approximation I(t) ~ M. 

(v) To carry out the average, run the steps 1-4 independently, say, v times, with 
v being a sufficiently large number, and take the arithmetic mean. 

Note that there is no need to recalculate the value r m i n after each tunneling: 
this should be done only if the pair with r min reacts. 

3.2. Recombination in the presence of diffusion 

In the general case when the atoms Hp and He not only recombine, but also diffuse, 
the algorithm becomes more sophisticated. Usually, we would consider diffusion to 
occur by microscopic random jumps according to the law dl = lo y/ D dt, where dl 
is the length of the random jump, lo is a random (isotropic) direction, and D is 
the diffusion coefficient. However, since the time between individual recombination 
events may be very large compared to dt, a huge number of diffusion jumps would 
be required to simulate the recombination dynamics. 

To accelerate this algorithm, we have suggested in [14] a Random Walk on 
Spheres based algorithm (more details about the Random Walk on Spheres method 
can be found in Ref. rat The idea behind this method is simple. Around each 
diffusing atom, we construct a disk of maximal radius which does not contain any 
atom Hp or He . We then simulate the random exit times T/ ; of the atom from 
these disks. We chose the atom which has a minimal exit time, and let this atom 
jump out of the disk, so that the new random position of the atom is uniformly 
distributed on the boundary of this disk. The distribution of the exit time is known 
(e.g., see [14], and we are thus able to simulate the random time dt according to 
this distribution, giving us the time t t + dt. 

This method has shown high efficiency for calculation of the intensity, as 
reported in P30- For our purpose however this approach cannot be directly applied: 
for viewing at the aggregation process where clusters of Hp atoms are separated 
by the clusters of H c atoms, we need to simulate the process step by step, with 
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sufficiently small time intervals. Therefore, we have introduced a mesh, and the 
atoms were diffusing over the mesh with the frequency A 5 = Dn for Hp , and 
A 6 = Dp for He • 

The general code described above remains the same. We have only to sample 
an additional event that the atom Hp makes a jump with frequency A 5 , and the 
atom He makes a jump with frequency Ag. Then, taking A 5 and Ae as described 
above, put A = ^ j =1 A*, and calculate the probabilities as p* = Aj/A, i — 1 ,... 6 . 

Analogously the resorption and exchange of He and Hp positions are 
simulated, using the rates Wh p , Wh c , and ape, olcpi respectively, which introduces 
four more events. Finally, the event of appearance a new Hp atom is simulated 
according to the accretion rate, or flux F, of H atoms onto the surface of a dust 
grain. Thus in total we have 11 events to be sampled at each time step. 

4. Simulation results 

As mentioned in the introduction, inhomogeneous fluctuations lead to the 
formation of clusters. The clustering slows down the reaction considerably, because 
only particles near the boundary of the clusters are likely to react, while particles 
inside the cluster have to diffuse to the boundary before they have a chance to 
react with a particle of the other type. In other words, fluctuations induce the 
formation of a mosaic of continuously growing domains which contain only one of 
the two components, Hp or He . 

Our simulations are adapted to take place on a spatial and temporal scale 
comparable to that for the surface of interstellar dust grains, given in [5j and [ 2 j. 
In Figure 1 we show 6 samples of random concentration fields of Hp and He when 
the system reaches a stochastically stationary regime. The segregation phenomen 
is clearly seen: large clusters of Hp and He are negatively correlated, and should 
be taken into account when evaluating the concentration of the molecular hydrogen 
from the equation Q. Without this correction, the model would much overestimate 
the efficiency rj given by (J5]) . The evaluation of ?/ in the framework of this model 
will be done in the forthcoming paper. 
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Figure 1 . Six samples of Hp and He atoms during the stochastically 
stationary regime on a grain of radius 0.1 gm. It is seen, that a segregation 
is formed with a characteristic distance between large clusters of Hp and 
He atoms. 
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5. Summary and conclusion 

Stochastic model and Monte Carlo simulation algorithm are constructed for solving 
a nonlinear system of inhomogeneous 2D Smuluchowski equations with random 
source term for simulation of H 2 formation on grain surfaces. The general system 
of inhomogeneous Smoluchowski type equations is used to govern the formation 
of H 2 from two hydrogen atoms on the surface of interstellar dust grains with 
physisorption and chemisorption sites. Both tunneling and diffusion mechanisms 
are taken into account. We focus in this study on the spatial distribution, and 
numerically investigate the segregation in the case of a source with a continuous 
generation in time and randomly distributed in space. The stochastic particle 
method presented is based on a probabilistic interpretation of the underlying 
process as a stochastic Markov process of interacting particle system in discrete 
but randomly progressed time instances. The segregation is analyzed through the 
correlation analysis of the vector random held of concentrations which appears to 
be isotropic in space and stationary in time. Note that the suggested model can be 
easily extended to more general situations, in particular, additional capture centers 
may be introduced in the system where the Hp and He may recombine without 
giving a contribution to the molecular hydrogen. This model will be presented in 
the forthcoming paper. 
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